function var_out  = parFileDSS(simType, avgCorr, varCorr, usebihist, inputPath, varIn, noSim, outputFilePath, krigType,...
    secVar, bihistFile, bounds, XX, YY, ZZ, localCorr, globalCorr, auxVar, krigANG, krigRANGE,...
    varANG,varRANGE, varType, varNugget,zoneFileName,noZones,noClasses,rescale,lvmFile, usePseudoHard, ...
    pseudoHardFile, pseudoCorr, covTab)
%% Run external DSS exectuable and return realization
%
% [input]
%   simType (string) - SGEMS or GEOEAS type
%   avgCorr (0/1) - correct mean
%   varCorr (0/1) - correct variance 
%   usebihist (0/1) - use joint probability distributions
%   inputPath (string) - path for folder with input data and DSS executable
%   varIn (string) - name of the file with input experimental data
%   noSim - number of realizations, 
%   outputFilePath (string) - full path for output file name without extension
%   krigType (0/1/2/3/4/5) - kriging type
%   secVar (string) - fullpath for secondary variable
%   bihistFile (string) - fullpath for joint distribution file
%   bounds (noZones x 2) - min and max of variable to be simulated
%   XX (1 x 3) - Number of cells, origin, size in XX
%   YY (1 x 3) - Number of cells, origin, size in YY
%	ZZ (1 x 3) - Number of cells, origin, size in ZZ
%   localCorr (string) - fullpath for collocated correlation coefficient
%       file
%   globalCorr - global correlation coefficient between 2 variables
%   auxVar string) - fullpath for sec variable for joint distribution
%   krigANG (1 x 3) - azimuth, rake and dip
%   krigRANGE (1 x 3) - major, minor, vertical
%   varANG (1 x 3) - azimuth, rake and dip
%   varRANGE (noZones x 3) - major, minor, vertical
%   varType(noZones x 1) - variogram type
%   varNugget (noZones x 1) - nugget effect
%   zoneFileName(string) - fullpath for zone file 
%   noZones - number of zones in zonefile
%   noClasses - number of classes for joint distribution
%   rescale (0/1) - rescale secondary variable for co-simulation
%   lvmFile (string) - fullpath for local varying means
%   usePseudoHard (0/1) - for point distributions
%   pseudoHardFile (string) - fullpath for local pdfs
%   pseudoCorr 0/1) - correct local pdfs
%   covTab (1 x 3) - size of the covariance matrix to be stored in mem
% 
% [output]
%   var_out - realization
%
%
% Leonardo Azevedo - 2019
% CERENA/Instituto Superior T�cnico (Portugal)
%

fid=fopen([inputPath,'ssdir.par'],'w');
fprintf(fid, '%s \n','#*************************************************************************************#');
fprintf(fid, '%s \n','#                                                                                     #');
fprintf(fid, '%s \n','#             PARALLEL DIRECT SEQUENCIAL SIMULATION PARAMETER FILE                    #');
fprintf(fid, '%s \n','#                                                                                     #');
fprintf(fid, '%s \n','#*************************************************************************************#');
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','#         Remember, # - Comment; [GROUP] - parameter group; CAPS - parameter          #');
fprintf(fid, '%s \n','#         Also, no space allowed in paths/filenames                                   #');
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','# here we define the hardata parameters');
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '\n');
fprintf(fid, '%s \n','[ZONES]');
fprintf(fid, '%s \n',['ZONESFILE = ',zoneFileName,'         # File with zones']);    
fprintf(fid, '%s \n',['NZONES = ', mat2str(noZones),'    # Number of zones']);
fprintf(fid, '\n');      
for n=1:noZones
    fprintf(fid, '%s \n',['[HARDDATA',num2str(n),']']);
    fprintf(fid, '%s \n',['DATAFILE  =',varIn(n,:),'   # Hard Data file']);
    fprintf(fid, '%s \n','COLUMNS   = 4                 # Number of columns in the data file (to eliminate)');
    fprintf(fid, '%s \n','XCOLUMN   = 1                 # Column number of X coordinate');
    fprintf(fid, '%s \n','YCOLUMN   = 2                 # Column number of Y coordinate');
    fprintf(fid, '%s \n','ZCOLUMN   = 3                 # Column number of Z coordinate');
    fprintf(fid, '%s \n','VARCOLUMN = 4                 # Column number for the variable');
    fprintf(fid, '%s \n','WTCOLUMN  = 0                 # Column number for the declustering weight');
    fprintf(fid, '%s \n',['MINVAL   = ', mat2str(bounds(n,1)),'          # Minimun threshold value']);
    fprintf(fid, '%s \n',['MAXVAL   = ', mat2str(bounds(n,2)),'          # Maximum threshold value']);
    fprintf(fid, '%s \n','USETRANS  = 1                 # Check if transformation file will be used');
    fprintf(fid, '%s \n','TRANSFILE = Cluster.trn       # Transformation file');
end       
fprintf(fid, '\n');      
fprintf(fid, '%s \n','[HARDDATA]');
fprintf(fid, '%s \n',['ZMIN     = ',  mat2str(min(bounds(:,1))),'            # Minimum allowable data value ']);
fprintf(fid, '%s \n',['ZMAX     = ',  mat2str(max(bounds(:,2))),'            # Maximum allowable data value']);
fprintf(fid, '%s \n','LTAIL     = 1                   # Specify the back-transform implementation in the lower tail ');
fprintf(fid, '%s \n',['LTPAR    = ',  mat2str(min(bounds(:,1))),'            # Parameter for the ltail=2']);
fprintf(fid, '%s \n','UTAIL     = 1                   # Specify the back-transform implementation in the upper tail ');
fprintf(fid, '%s \n',['UTPAR    = ',  mat2str(max(bounds(:,2))),'            # Parameter for the utail=2']);
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','# here we define some parameters for the simulation');
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','[SIMULATION]');
fprintf(fid, '%s \n',['OUTFILE   = ', outputFilePath,'         # Filename of the resultant simulations']);
fprintf(fid, '%s \n',['NSIMS     = ',num2str(noSim),'   # Number of Simulations to generate']);
fprintf(fid, '%s \n','NTRY       = 10                   # Number of simulation for bias correction ');
fprintf(fid, '%s \n',['AVGCORR   = ', num2str(avgCorr), '       # Variance correction flag ']);
fprintf(fid, '%s \n',['VARCORR   = ', num2str(varCorr), '       # Average correction flag']);
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','# here we define the output grid (and secondary info grid)');
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','[GRID]');
fprintf(fid, '%s \n','# NX, NY and NZ are the number of blocks per direction ');
fprintf(fid, '%s \n',['NX        = ', num2str(XX(1,1))]);
fprintf(fid, '%s \n',['NY        = ',num2str(YY(1,1))]);
fprintf(fid, '%s \n',['NZ        = ',num2str(ZZ(1,1))]);
fprintf(fid, '%s \n','# ORIGX, ORIGY and ORIGZ are the start coordinate for each direction');
fprintf(fid, '%s \n',['ORIGX     = ', num2str(XX(1,2))]);
fprintf(fid, '%s \n',['ORIGY     = ', num2str(YY(1,2))]);
fprintf(fid, '%s \n',['ORIGZ     = ', num2str(ZZ(1,2))]);
fprintf(fid, '%s \n','# SIZEX, SIZEY and SIZEZ is the size of blocks in each direction ');
fprintf(fid, '%s \n',['SIZEX     = ', num2str(XX(1,3))]);
fprintf(fid, '%s \n',['SIZEY     = ', num2str(YY(1,3))]);
fprintf(fid, '%s \n',['SIZEZ     = ', num2str(ZZ(1,3))]);
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','# here we define some general parameters');
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','[GENERAL]');
fprintf(fid, '%s \n','NULLVAL       = -99999.25                                     # Definition of the null value');
fprintf(fid, '%s \n',['SEED         = ', num2str(int32(rand*1000000)),'             # Seed for the pseudorandom number generator - if you want repeatable simulations']);
fprintf(fid, '%s \n','USEHEADERS    = 1                                             # Seed for the pseudorandom number generator - if you want repeatable simulations');
fprintf(fid, '%s \n',['FILETYPE     = ',simType, '                                  # accepts GEOEAS and SGEMS. default is sgems file type. used for output ']);
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','# here we define the parameters for search');
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','[SEARCH]');
fprintf(fid, '%s \n','NDMIN   = 1                   # Minimum number of original data that should be used to simulate a node');
fprintf(fid, '%s \n','NDMAX   = 32                  # Maximum number of original data that should be used to simulate a node');
fprintf(fid, '%s \n','NODMAX  = 12                  # Maximum number of previouly simulated nodes to be used to simulate a new node');
fprintf(fid, '%s \n','SSTRAT  = 1                   # Two-part search / data nodes flag ');
fprintf(fid, '%s \n','MULTS   = 0                   # Multiple grid simulation flag ');
fprintf(fid, '%s \n','NMULTS  = 1                   # Number of multiple grid refinements ');
fprintf(fid, '%s \n','NOCT    = 0                   # Number of original data to use per octant');
fprintf(fid, '%s \n',['RADIUS1 = ', num2str(krigRANGE(1,1)),'                # Search radii in the maximum horizontal direction']);
fprintf(fid, '%s \n',['RADIUS2 = ', num2str(krigRANGE(1,2)),'                # Search radii in the minimum horizontal direction']);
fprintf(fid, '%s \n',['RADIUS3 = ', num2str(krigRANGE(1,3)),'                # Search radii in the vertical direction']);
fprintf(fid, '%s \n',['SANG1   = ', num2str(krigANG(1,1)),'                   # Orientation angle parameter of direction I (degrees)']); 
fprintf(fid, '%s \n',['SANG2   = ', num2str(krigANG(1,2)),'                   # Orientation angle parameter of direction II (degrees)']);
fprintf(fid, '%s \n',['SANG3   = ', num2str(krigANG(1,3)),'                   # Orientation angle parameter of direction III (degrees)']);
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','# here we define the kriging information, and secondary info when applicable');
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','[KRIGING]');
fprintf(fid, '%s \n',['KTYPE        = ', num2str(krigType),'                # Kriging type ;0=simple,1=ordinary,2=simple with locally varying mean,']);
fprintf(fid, '%s \n','# 3=external drif, 4=collo-cokrig global CC,5=local CC (KTYPE)');
fprintf(fid, '%s \n',['COLOCORR     = ', num2str(globalCorr),'                # Global CC to ktype=4 ']);
fprintf(fid, '%s \n',['SOFTFILE     = ', secVar,'                                  # Filename of the soft data']);
fprintf(fid, '%s \n',['LVMFILE      = ', lvmFile,'                                 # FOR KTYPE=2']);
fprintf(fid, '%s \n','NVARIL        = 1                                             # Number of columns in the secundary data file');
fprintf(fid, '%s \n','ICOLLVM       = 1                                             # Column number of secundary variable ');
fprintf(fid, '%s \n',['CCFILE       = ', localCorr,'                               # Filename of correlation file for local correlations (ktype=5)']);
fprintf(fid, '%s \n',['RESCALE      = ', num2str(rescale),'                         # rescale secondary variable (for ktype>=4)']);
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','# here we define the variogram to use. if more than 1, use [VARIOGRAM2]');
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
for n=1:noZones
    fprintf(fid, '%s \n',['[VARIOGRAMZ',num2str(n),']']);
    fprintf(fid, '%s \n',['NSTRUCT  = ',num2str(size(varRANGE,2)/4),'                # Number of semivariograms structures (NST(1))']);
    fprintf(fid, '%s \n',['NUGGET    =', num2str(varNugget(n,1)),'                                               # Nugget constant (C0(1))']);
    for i=1:size(varRANGE,2)/4
        fprintf(fid, '%s \n',['[VARIOGRAMZ',num2str(n),'S',num2str(i),']']);
        fprintf(fid, '%s \n',['TYPE =', num2str(varType(n,1)), '                   # Struture type ;1=spherical,2=exponential,3=gaussian (IT(i))']);
        fprintf(fid, '%s \n',['COV  =', num2str(varRANGE(n,4*i))  ,'               # C parameter "COV + NUGGET = 1.0" (CC(i))']);
        fprintf(fid, '%s \n',['ANG1 =', num2str(varANG(1,1))  ,'                   # Geometric anisotropy angle I (ANG1(i))']);
        fprintf(fid, '%s \n',['ANG2 =', num2str(varANG(1,2))  ,'                   # Geometric anisotropy angle II (ANG2(i))']);
        fprintf(fid, '%s \n',['ANG3 =', num2str(varANG(1,3))  ,'                   # Geometric anisotropy angle III (ANG3(i))']);
        fprintf(fid, '%s \n',['AA   =', num2str(varRANGE(n,4*i-3)),'               # Maximum horizontal range (AA(i))']);
        fprintf(fid, '%s \n',['AA1  =', num2str(varRANGE(n,4*i-2)),'               # Minimum horizontal range (AA1) ']);
        fprintf(fid, '%s \n',['AA2  =', num2str(varRANGE(n,4*i-1)),'               # Vertical range (AA2)']);
    end
end
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','# here we define parameters for joint DSS ');
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
for n=1:noZones
    fprintf(fid, '%s \n',['[BIHIST',num2str(n),']']);
    fprintf(fid, '%s \n',['USEBIHIST        = ', num2str(usebihist(1,1)),' #Use Bihist? 1-yes 0-no']);
    fprintf(fid, '%s \n',['BIHISTFILE       = ', bihistFile(n,:),'   # bihistogram file']);
    fprintf(fid, '%s \n',['NCLASSES         = ', num2str(noClasses),'                           # number of classes to use']);
    fprintf(fid, '%s \n',['AUXILIARYFILE    = ', auxVar,'                   # auxiliary image']);            
end
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','# here we define the debug parameters - probably you wont need this');
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','[DEBUG]  ');
fprintf(fid, '%s \n','DBGLEVEL  = 1                 # 1 to 3, use higher than 1 only if REALLY needed');
fprintf(fid, '%s \n','DBGFILE   = debug.dbg         # File to write debug');
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','# here we define parameters for COVARIANCE TABLE - reduce if memory is a problem ');
fprintf(fid, '%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid, '%s \n','[COVTAB]');
fprintf(fid, '%s \n',['MAXCTX = ', num2str(covTab(1,1))]);
fprintf(fid, '%s \n',['MAXCTY = ', num2str(covTab(1,2))]);
fprintf(fid, '%s \n',['MAXCTZ = ', num2str(covTab(1,3))]);
fprintf(fid,'%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid,'%s \n','# here we define parameters for BLOCK KRIGING - if you are not block kriging useblocks should be 0 ');
fprintf(fid,'%s \n','#-------------------------------------------------------------------------------------#');
fprintf(fid,'%s \n','[BLOCKS]');
fprintf(fid,'%s \n','USEBLOCKS  = 0   # 1 use, 0 no');
fprintf(fid,'%s \n','BLOCKSFILE = no file');
fprintf(fid,'%s \n','MAXBLOCKS  = 100');
fprintf(fid, '%s \n','[PSEUDOHARD]');
fprintf(fid,'%s \n',['USEPSEUDO  = ', num2str(usePseudoHard), '  # 1 use, 0 no  pseudo hard data is point distributions that are simulated before all other nodes']);
fprintf(fid,'%s \n',['PSEUDOFILE = ', pseudoHardFile,'           # file']);
fprintf(fid,'%s \n',['PSEUDOCORR = ', num2str(pseudoCorr),'               # correct simulated value with point']);
fclose(fid);
%%   RUN DSS
temp = [inputPath,'DSS.C.64.exe   ',inputPath,'ssdir.par'];

% hack to silent mode in command window
[T,Status,results] = evalc('system(temp)');

%%  IMPORT SIMULATIONS
if (strcmp(simType ,'SGEMS') == 1);
    
    varStrut    = sgems_read([outputFilePath,'_1.sgems']);
    var_out     = varStrut.data;
    clear varStrut
    
else
    [~,var_out]=import_gslib([outputFilePath,'_1.out']);
end

end